% this is the script for plotting MRA data
% Version: 2.2
clc;
clear all;
close all;

%Read data file
    %fileName=input('Input a file name: ','s');
    fileName='test_grain_list.txt';
    f_id=fopen(fileName);
    label=textscan(f_id,'%s %s %s %s %s %s %s %s',1);
    data=fscanf(f_id,'%f',[8 inf]);
    f_return=fclose(f_id);
    
    id=data(1,:)';
    area=data(2,:)';
    density=data(3,:)';
    sizeDens=data(4,:)';
    clusterDens=data(5,:)';
    freq=data(6,:)';
    de=data(7,:)';
    nFreq=data(8,:)';
    
    diameter=sqrt(area/pi)*2;
    diameter=diameter';
    %minDiameter=min(diameter);
    logDiameter=log10(diameter);

%extract data

    
    numSizeClass=3;
    sizeClass(1:size(id,1))=0;
    sizeClass=sizeClass';
    sizeInt=0.0000001+(max(logDiameter)-min(logDiameter))/numSizeClass;

    gridIntv=0:0.005:1;
    [coX coY]=meshgrid(gridIntv,gridIntv);
    
    for k=1:numSizeClass
        selected=find(logDiameter>=(sizeInt*(k-1)+min(logDiameter)) & logDiameter<(sizeInt*k+min(logDiameter)));
        sizeClass(selected)=k;
        numClassParticles(:,k)=size(selected,2);
        x=(nFreq(selected));
        y=(density(selected));
        z=(logDiameter(selected));

        fitLines(:,k)=polyfit(x,y,1);

        for i=1:size(gridIntv,2)
            for j=1:size(gridIntv,2)
                coZ(i,j,k)=size(find(power(x-coX(i,j),2)+power(y-coY(i,j),2)<power(0.01,2)),1);
                coZ(i,j,k)=coZ(i,j,k)/size(id,2);
            end;
        end;
    end;

    x=nFreq;
    y=density;
    z=logDiameter;

    k=k+1
    numClassParticles(:,k)=size(id,1)
    fitLines(:,k)=polyfit(x,y,1)

    for i=1:size(gridIntv,2)
        for j=1:size(gridIntv,2)
            coZ(i,j,k)=size(find(power(x-coX(i,j),2)+power(y-coY(i,j),2)<power(0.01,2)),1);
            coZ(i,j,k)=coZ(i,j,k);
        end;
    end;

%plot data
    %hold on;
    
    %intv=0:0.001:1;
    %rangeMaxX=max(max(x), max(y))+0.1;
    %rangeMinX=min(min(x), min(y))-0.1;
    %rangeMaxY=max(max(x), max(y))+0.1;
    %rangeMinY=min(min(x), min(y))-0.1;
    rangeMaxX=max(nFreq)+0.05;
    rangeMinX=min(nFreq)-0.05;
    rangeMaxY=max(density)+0.05;
    rangeMinY=min(density)-0.05;
    
    intvX=rangeMinX:0.001:rangeMaxX;
    intvY=rangeMinY:0.001:rangeMaxY;
    [XI,YI] = meshgrid(intvX, intvY);
    
    axis_font_size=13;
    conLines=10;
    plotNum=[2 3];
    
   
    subplot(plotNum(1),plotNum(2),1);
    ZI = griddata(coX,coY,coZ(:,:,4),XI,YI);
    contour(XI,YI,ZI,conLines); 
    grid on;
    axis equal;
    axis([rangeMinX rangeMaxX rangeMinY rangeMaxY]);
    xlabel('Neighboring Density','fontsize',axis_font_size);
    ylabel('Area Density','fontsize',axis_font_size);
    title(['Total (',num2str(numClassParticles(4)),')'],'fontsize',16);
    colorbar;
    
    subplot(plotNum(1),plotNum(2),[2 3]);
    %ZI = griddata(coX,coY,coZ(:,:,4),XI,YI);
    %mesh(XI,YI,ZI); 
    scatter3(x,y,z,diameter*5,sizeClass,'filled');
    %axis([rangeMinX rangeMaxX rangeMinY rangeMaxY]);
    xlabel('Neighboring Density','fontsize',axis_font_size);
    ylabel('Area Density','fontsize',axis_font_size);
    zlabel('Normalized Size (log)','fontsize',axis_font_size);
    %colorbar;
    view(-50,45);


    
    subplot(plotNum(1),plotNum(2),4);
    ZI = griddata(coX,coY,coZ(:,:,1),XI,YI);
    contour(XI,YI,ZI,conLines); 
    grid on;
    axis equal;
    axis([rangeMinX rangeMaxX rangeMinY rangeMaxY]);
    title(['Small (',num2str(numClassParticles(1)),')'],'fontsize',16);
    colorbar;

    subplot(plotNum(1),plotNum(2),5);
    ZI = griddata(coX,coY,coZ(:,:,2),XI,YI);
    contour(XI,YI,ZI,conLines); 
    grid on;
    axis equal;
    axis([rangeMinX rangeMaxX rangeMinY rangeMaxY]);
    title(['Medium  (',num2str(numClassParticles(2)),')'],'fontsize',16);
    colorbar;

    subplot(plotNum(1),plotNum(2),6);
    ZI = griddata(coX,coY,coZ(:,:,3),XI,YI);
    contour(XI,YI,ZI,conLines); 
    grid on;
    axis equal;
    axis([rangeMinX rangeMaxX rangeMinY rangeMaxY]);
    title(['Large  (',num2str(numClassParticles(3)),')'],'fontsize',16);
    colorbar;
    
    %scatter3(x,y,z,'filled');
    %scatter(x,y,diameter(selected)*3,diameter(selected),'filled');
    %scatter(nFreq,density);
    %scatter(nFreq,logDiameter,area);

    %hold off;    

    %colormap jet;
    
    %axis ([0 1 0 1]);
    %axis ([min(x) max(x) min(y) max(y)]);
    %zlabel('Normalized Diameter (log)','fontsizeDensDens',axis_font_sizeDensDens);
    %grid off;
    %camlight left;
